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Abstract 


The present study addresses an exponentially fitted finite difference method 
to obtain the solution of singularly perturbed two-point boundary value 
problems (BVPs) having a boundary layer at one end (left or right) point 
on uniform mesh. A fitting factor is introduced in the derived scheme 
using the theory of singular perturbations. Thomas algorithm is employed 
to solve the resulting tri-diagonal system of equations. The convergence of 
the presented method is investigated. Several model example problems are 
solved using the proposed method. The results are presented with terms of 
maximum absolute errors, which demonstrate the accuracy and efficiency 
of the method. It is observed that the proposed method is capable of 
producing highly accurate results with minimal computational effort for a 
fixed value of step size h, when the perturbation parameter tends to zero. 
From the graphs, we also observed that a numerical solution approximates 
the exact solution very well in the boundary layers for smaller value of «. 
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1 Introduction 


Singular perturbation problems are of mainly deal in fluid mechanics and 
other areas of practical/applied mathematics. The solution of the singularly 
perturbed boundary value problems (BVPs) has a multi-scale nature. The 
solution varies rapidly in some parts of the domain and varies slowly in some 
other parts of the domain. The numerical solution of singular perturbation 
problems (SPPs) is far from trivial, because of the boundary layer behav- 
ior of the solution. There are many physical situations in which the sharp 
changes occur inside the domain of interest, and the narrow regions across 
which these changes take place are usually referred as Navier-Stokes flow 
problems, involving high Reynolds number [4, 17, 28], mathematical models 
of liquid crystal materials and chemical reactions, control theory, and elec- 
trical networks [6, 7, 30]. These quick shifts can be managed by fast scales, 
magnified scales, or stretched scales, but not by slow scales. The domain of 
integration is typically divided into two subdomains, with a distinct scheme 
being applied to each subdomain as a common approach to solving this type 
of problem. In recent years, a large number of analytical methods have 
been proposed (see [22, 21, 2, 20, 19, 11, 16, 29]. Numerical methods based 
schemes with and without fitting factors on boundary value techniques and 
initial value techniques are given in [9, 1, 12, 13, 23, 14]. Phaneendra and 
Lalu[24] presented Gaussian quadrature for two-point singularly perturbed 
BVPs with the exponential fitting with a layer at one endpoint, dual bound- 
ary layers, and internal boundary layers. In this paper, the given BVP is 
reduced into an equivalent first-order differential equation with the pertur- 
bation parameter as a deviating argument. Then, the Gaussian two-point 
quadrature technique with exponential fitting is implemented to solve the 
first-order equation with deviating parameters. Mishra and Saini [18] stud- 
ied the Liouville-Green transform to solve a singularly perturbed two-point 
BVP with a right-end boundary layer. Articles [3, 5, 8, 9, 31] proposed differ- 
ent numerical approaches combining fitted mesh methods and fitted operator 
methods employed by several researchers for solving SPPs, whereas Kadal- 
bajoo and Kumar [10] presented a detailed outline on the numerical methods 
for solving SPPs. Indeed these existing numerical methods are mostly based 
on fitted operator techniques or use reasonable theoretical information re- 
garding the solutions, which forms a limitation of these approaches. An 
efficient method of numerical integration for a class of singularly perturbed 
two-point BVPs at one endpoint (either left or right) has been presented in 
[25]. Ranjan, Prasad, and Alam [27] developed a simple method of numerical 
integration for a class of singularly perturbed two-point BVPs at one end- 
point (either left or right). Ranjan and Prasad [26] proposed a fitted finite 
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difference scheme for solving singularly perturbed two-point BVPs having 
boundary layer at left or right endpoints. Madhu Latha, Phaneendra, and 
Reddy [15] presented a numerical solution of SPP using numerical integration 
with an exponential fitting factor. 


In view of the wealth of literature on SPPs, we raise the question of 
whether there are other ways to attack SPPs, namely ways that are very 
easy to use and ready for computer implementation. In response to this 
need for a fresh approach to SPPs, we propose and illustrate in this paper 
a fitted finite difference technique for singularly perturbed two-point BVPs 
with a boundary layer on the left (or right) end of the underlying interval. 
Numerical experience with several linear examples is described. 


The paper is organized as follows: Section 2 presents the description of the 
presented new effective method to solve a second-order singularly perturbed 
two-point BVP. In Section 3, the convergence of the presented method is 
investigated. To demonstrate the accuracy and efficiency of the presented 
method, numerical experiments are carried out for several model test prob- 
lems, and the results are shown in tables in Section 4. Finally, the discussions 
and conclusions are presented in the last section 5. 


2 Statement of the problems 


Consider the singularly perturbed two-point BVPs of the following type: 
ev" (t) + r(t)u'(t) + s(t)v(t) = W(t) on N = (0, 1], (1) 
subject to the boundary conditions and interval conditions, 
v(0) =a, v(1) =, (2) 


where ¢ is a small positive perturbation parameter (0 < ¢ << 1). Further- 
more, the functions r(t), s(t), and w(t) are continuously differentiable func- 
tions in [0,1], where a and @ are constant. In the scenario where we assume 
that r(t) > M > 0 holds true for the entire interval [0,1], with M repre- 
senting a positive constant, the boundary layer is expected to occur in the 
vicinity of t = 0. On the other hand, if we consider that r(t) <M < 0 holds 
throughout the interval [0,1], with M being a negative constant, then the 
boundary layer is anticipated to be located near t = 1. 
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2.1 Description of the method for left-end boundary 
layer problems 


In this subsection, we describe the proposed method for the solution of the 
problem (1)—(2) having boundary layer at left-end point of the interval con- 
sidered. 


The solution of (1) with (2) is of the following form (see . 22—261[22]): 
0 _— rt(r@_ st 
v(t) = wolt) + 2 (ao — unl BPFH 4 of), 8) 


where v(t) denotes the simplified problem’s solution: 
r(t)up(t) + s(t)vo(t) =v), —-vo(1) = B. (4) 


By considering the Taylor series expansions of r(t) and s(t) around the point 
t = 0 up to their respective first terms, we can simplify (3) as follows: 


u(t) = v9(t) + (a — v0(0)) ee FO)* + o(e). (5) 


Taking the limit as h — 0 and applying (3) to the point t = t; = th,i = 
0,1,2,...,N, we obtain 


r2(0)—es(0) 


lim v(th) = v9 (0) + (ao — v0(0)) a si )ie + o(€), (6) 


h-0 


where p = h/e, the first and second-order approximations have been used as 
below: 


ge 3Vj41 — 20; — Vi-1 


7 
u Ah ? ( ) 
U5 = 2u; + Uj 
an an (8) 
Substituting (7) and (8) in (1), we have 
Vig-1 — 205 + Vi-1 3Uj41 — 20; — Vi-1 
E | 2 | +7; | Th + 5,0; = V4. (9) 


Introducing the fitting factor o(p) into the aforementioned approach, we 
obtain the following result: 


7 


as — 2u; + “=| i 


3Ui41 = 20; = 
h2 


Vi-1 = 
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The determination of the fitting factor o(p) aims to ensure that the solution of 
the difference scheme described in (10) achieves uniform convergence towards 
the solution of (1) with (2). 


By multiplying (10) by # and considering the limit as h —> 0, the result 
of (10) is as follows: 


0 
2 [usta 2U; U;-1] t 0) [Bui41 2U; U;-1] = 0. (11) 


Let p = Co) s0(0), By using (6), we get 


limpo (v(ih — h) + v(ih + h) — 2v(ih)) = (ao — v9(0)) e~h"? (eM? + eH? — 2), 
limp_o (3u(ih + h) — 2v(éh) — v(ih — h)) = (ao — vo(0)) eM"? (Be7 HK”? — 2 — eM?) . 


By using the above equations in (11), we get 


“we a (‘ os) ") aoe. (12) 


which is a required fitting factor o (p) . 
Finally, from (11) with the value of o(p) given by (12), we obtain the following 
three-term recurrence relationship: 


Py Oe, ty Sy (1,2, 8, xan = 1), (13) 
where 
es 
Qi = FF + ah ~ Sis 


Ay = yi. 


Equation (13) generates an (N — 1) equations system involving (N — 1) un- 
determined ranging from v; to vy_;. These (N—1) equations together with 
the boundary conditions equation (2), are sufficient to solve the obtained 
tri-diagonal system with the help of an efficient solver called the Thomas al- 
gorithm, commonly called as the “Discrete Invariant Imbedding algorithm”. 


2.2 Description of the method for right-end boundary 
layer problems 


In this subsection, we will describe the proposed method for the solution of 
the problem (1)—(2) having boundary layer at right-end point of the interval 
considered. 
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The solution of (1) with (2) is of the following form . 22—261[22]): 
v(t) = wolt) + A (ao = wp(ayyen BEP“HBM 4 of), (48 
where yo(t) denotes the simplified problem’s solution: 
r(t)up(t) + s(t)vo(t) = v(t), —-vo(1) = B. (15) 


By considering the Taylor series expansions of r(t) and s(t) around the point 
t = 0 up to their respective first terms, we can simplify (14) as follows: 


u(t) = vo(t) + (ao — v9(0)) e7® ~ FaI)* + fe). (16) 


Taking the limit as h — 0 and applying (3) to the point t = t; = th,i = 
0,1,2,...,N, we obtain 


limp_so v(ih) = v9 (0) + (ao —= yo(0)) e c Das Jee + o(e), (17) 


where p = h/e. 
After multiplying (10) by h and taking the limit as h > 0, (10) converts 
into the following form: 


0 
: [usta — 2u; - Ui-1| + m ) [Bui441 20; U;-1] = 0. (18) 


r?(0)—es(0) 


Let w= 7(0) 


. By using (17), we get 


limpso (v(ih — h) + v(ih + h) — 2v(ih)) = (ao — v9(1)) eM? (eM? + eH? — 2), 
limp_s0 (3u(éh + h) — 2u(ih) — v(ih — h)) = (ap — v9(1)) eM"? (3e74? — 2 — eM?) . 


By substituting the aforementioned equations into (18), we get 


ei a oa (ee = es(1)) "| r(0)p (19) 


2r(1) 4? 
which is a required fitting factor o (~) for right-end boundary layer problem. 
3 Convergence analysis 


This section focuses on the analysis of the convergence of the method. 


Definition 1 (Consistency). Let 


od; |v] = Lpnv(ti) — Lev(ti), 4 Dds saa lV. 
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In this context, v denotes a smooth function defined on the interval J = [0, 1], 
and L,, represents the discrete difference operator. Consequently, the differ- 
ence equation (13)—(2) exhibits consistency with the corresponding differen- 
tial equation (1)—(2), if 


|pi[v]] > 0Oash—-o. 


The quantities ¢;[v],¢ = 1,2,3...,N is called the local truncation (or local 
discretization) errors. 


Definition 2. The differential equation (13)—(2) is said to possess local pth- 
order accuracy when, for suitably smooth data, there exists a positive con- 
stant C that remains independent of h and « such that 


maxi<i<n|¢i[v]| < Ch?. 


The agreement between the differential equation (13)—(2) and (1)-(2), 
along with its locally second-order accuracy, is established through the lemma 
provided below. 


Lemma 1. If v € C?(Z), then 


ih 
ld; [v]] < maxes, .<t<te44 {ene 10h), (5 1;2,3.0<,N = 1: 


Proof. By definition, 


- Vig — 205 + Vi-1 v 3Uj41 — 2U; — Vi-1 
= v; t : 
h2 : Ah 
h?.. ht. h h? 
i= a2 | Bol + out +h | nf el | mu seo} 
oeh? aa rjh VW 
lor — ee 12 [vi | - Neco ne 4 Iv , 
r;h 
ols amas {pti} +00?) 


ldi| << O(h). i1=1,2,3,...,N—1. 


As a result, the intended outcome is attained. 


We will now examine the proposed method’s convergence across the entire 
interval range 0 < t < 1. We write the tridiagonal system (13) in the matrix- 
vector form 

WV =D, (20) 


where W = (a;;), 1 < 1,7 < N-—1 is a tridiagonal matrix of order N — 1 
with 


Tran. J. Numer. Anal. Optim., Vol. 13, No. 4, 2023,pp 711-727 


Kumar, Kumar Sinha and Ranjan 718 


rjh 
Qj ji-1 = OE — , 
re 
Qj, = —20€6 — “4 §:h?, 
3hr; 
Aii41 = OE + T° 


and D = (d;) is a column vector with d; = h?¢,; for i = 1,2,3,..., N—1 with 
local truncation error ¢;: 


|oi| < O(A). (21) 


We also have 

WV — o(h) =D, (22) 
where V = (Vo, Vi, Va, V3,-.., Vw)’ and $(h) = (¢1(h), b2(h), o3(h),..., on (h))* 
stands for the local truncation error and the real solution, respectively. (20) 
and (22) give us 

W(V —V)= o(h). (23) 


Thus the error equation is 
WE = ¢(h), (24) 


where E = V —V = (eo, €1, €2,-.., en)’. If S¥ is the total of the components 
in the 7th row of W, then 


x N-1 _ -o¢e T1 
St = dja ay = 5 — Bt, 


* _yN-1 __ =agE 3rn-1 
Sn-1 = ast an-1,3 = 2 — gy + $N-1; 
Ss = j=l aj = Si = Bio. 


Since 0 < ¢ << 1, The matrix W is irreducible and monotone for sufficiently 
small h. As a result, W~! must exist and contain nonnegative elements. 
Therefore, we have from (24) that 


E=W~'6(h), (25) 
El] < Wo" NeA)I- (26) 


Let Gj; represent the (ki)th components of W~+. Since a; > 0, we have 
from the operations on matrices: 


N-1 


GyiS* =1, k=1,2,...,N—1. (27) 


j=l 
Therefore, its follows 


sil 1 1 1 


s Aki < zi = = ’ 
1 ming<i<w-159; Bio ~ |Biol 
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for some ig between 1 and N — 1, and Big = qG. 
Therefore, from (21), (25), and (27), we get 


N-1 
i=1 
which implies 
h 
ae “ 2, oN A. (30) 
Or 


Consequently, by applying the definitions and (29), we obtain: 
||| = Oh). 


This implies that the purposed method is the first-order rate of convergence 
on uniform mesh. As above, we can apply the same procedure for showing 
the purposed method is of first-order rate of convergence on uniform mesh 
for the right layer problem. 


4 Numerical illustrations 


The effectiveness of the purposed method has been demonstrated by imple- 
menting it on the three linear SPPs at left-end boundary layer as well as 
one problem involving a right-end boundary layer and presented the compu- 
tational results in the tables in terms of the maximum absolute errors EN. 
These examples have been chosen because they have been widely discussed in 
literature. For various values of mesh point N and perturbation parameter ¢, 
the EN are defined by EX = pax lleG) —v,|], where v(t;) and v; denote 
the exact and approximate solution, respectively. The double mesh principle 
is used to calculate the rate of convergence defined as r¥ = log, (er) . The 
purposed method is capable of achieving uniform results, when perturbation 
parameter ¢ tends to 0 for any fixed value of the mesh size h. 


Example 1. First, consider the following homogeneous SPP from [15]: 
ev" (t) + u' (t) — v(t) =0,  t € [0,1], 


with boundary condition v(0) = 1 and v(1) = 1. 
The exact solution is given by 


where m, = ————_——_——_ and mz = ———______. 


(-14+ V1 +4 4e) 
2€ 
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The maximum absolute errors for various values of N and singular per- 
turbation parameter € are presented in Table 1 for example 1. It can be easily 
observed from Table 1 that the maximum absolute errors tends uniformly, 
when the singular perturbation parameter ¢ tends to 0, for any fixed value of 
N = 1/h. Also, rates of convergence presented in Table 1 show that the pur- 
posed scheme is capable of producing almost first-order accurate uniformly 
convergent solution. In Figure 1, we present our solution and the exact solu- 
tion for various values of € and a fixed value of N. Clearly, as shown in the 
figure, the numerical solution and the exact solution are very close within the 
boundary layers for smaller values of ¢. 


Example 2. Consider the following non-homogeneous SPP involving a con- 
stant term f(t) [25, 15]: 


ev’ (t)+u'(t)=2, te [0,1], 


with boundary condition v(0) = 1 and v(1) = 1. The exact solution is given 
1-e!/® 


The maximum absolute errors for various values of N and singular per- 
turbation parameter € are presented in Table 2 for example 2. It can be easily 
observed from Table 2 that the maximum absolute errors tends uniformly, 
when the singular perturbation parameter ¢ tends to 0, for any fixed value 
of N = 1/h. In Figure 2, we present our solution and the exact solution for 
various values of ¢ and a fixed value of N. Clearly, as shown in the figure, the 
numerical solution and the exact solution are very close within the boundary 
layers for smaller values of e. 


Example 3. Consider the following non-homogeneous SPP involving a vari- 
able term f(t) (25, 15): 


ev" (t)+ u(t) =14 2t, te [0,1], 


with boundary condition v(0) = 1 and v(1) = 1. The exact solution is given 


by v(t) = 12895 (Qe — 1) + t( + 1 — 26). 


The maximum absolute errors for various values of N and singular per- 
turbation parameter ¢ are presented in Table 3 for example3. It can be easily 
observed from Table 3 that the maximum absolute errors tends uniformly, 
when the singular perturbation parameter ¢ tends to 0, for any fixed value of 
N =1/h. Also, rates of convergence presented in Table 3 show that the pur- 
posed scheme is capable of producing almost first-order accurate uniformly 
convergent solution. In Figure 3, we present our solution and the exact solu- 
tion for various values of € and a fixed value of N. Clearly, as shown in the 
figure, the numerical solution and the exact solution are very close within the 
boundary layers for smaller values of ¢. 
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Example 4. Lastly, consider the following homogeneous SPP at right-end 
boundary layer [25, 26, 15]: 


ev" (t) — v'(t) — (1+ e)v(t) =0, te [0,1], 


with boundary condition v(0) = 1+ eaxp(—(1+¢)/e) and v(1) = 1+1/e. The 
exact solution is given by v(t) = eA te)(@-D/e + et, 


The maximum absolute errors for various values of N and singular per- 
turbation parameter € are presented in Table 4 for Example 4. It can be 
easily observed from the Table 4 that the maximum absolute errors tends 
uniformly, when the singular perturbation parameter € tends to 0, for any 
fixed value of N = 1/h. Also, rates of convergence presented in Table 4 show 
that the purposed scheme is capable of producing almost first-order accurate 
uniformly convergent solution. In Figure 4, we present our solution and the 
exact solution for various values of ¢ and a fixed value of N. Clearly, as 
shown in the figure, the numerical solution and the exact solution are very 
close within the boundary layers for smaller values of ¢. 


Table 1: Computational results in terms of maximum absolute errors for various values 
of N and ¢ and the rate of convergence rN for Example 1 


N 16 32 64 128 256 512 
e€=10-> 0.0112 0.0057 0.0029 0.0014 0.0007 0.0003 
pal 0.9745 0.9749 1.0506 1.0000 1.2224 
e=10-° 0.0112 0.0057 0.0029 0.0014 0.0007 0.0003 
rN 0.9745 0.9749 1.0506 1.0000 1.2224 
e€=10-7 0.0112 0.0057 0.0029 0.0014 0.0007 0.0003 
rN 0.9745 0.9749 1.0506 1.0000 1.2224 
e=10-% 0.0112 0.0057 0.0029 0.0014 0.0007 0.0003 
rea 0.9745 0.9749 1.0506 1.0000 1.2224 
e=10-9 0.0112 0.0057 0.0029 0.0014 0.0007 0.0003 
ea 0.9745 0.9749 1.0506 1.0000 1.2224 


Table 2: Computational results in terms of maximum absolute errors for various values 
of N and ¢ for Example 2 


N 16 32 64 128 256 512 


é€=10-°> 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 
e=10-° 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 
e=10-" 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 
e=10-8 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 
e=10-? 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 
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Table 3: Computational results in terms of maximum absolute errors for various values 
of N and « and the rate of convergence r% for Example 3 


16 32 64 128 256 512 


0.0586 0.0303 0.0154 0.0079 0.0039 0.0019 
0.9516 0.9764 0.9630 1.0184 1.0375 
0.0586 0.0303 0.0154 0.0077 0.0039 0.0019 
0.9516 0.9764 1.0000 0.9814 1.0375 
0.0586 0.0303 0.0154 0.0078 0.0039 0.0019 
0.9516 0.9764 0.9814 1.0000 1.0375 
0.0586 0.0303 0.0154 0.0078 0.0039 0.0019 
0.9516 0.9764 0.9814 1.0000 1.0375 
0.0586 0.0303 0.0154 0.0078 0.0039 0.0019 
0.9516 0.9764 0.9814 1.0000 1.0375 


o a) oy o ia) 
II II lI lI lI 
3 3 3 3 
| 4 
oa i, a 
Ke} Go _ fon) oy 


OTe a 


Table 4: Computational results in terms of maximum absolute errors for various values 
of N and ¢ and the rate of convergence rN for Example 1 


16 32 64 128 256 512 


0.0112 0.0057 0.0029 0.0014 0.0007 0.0003 
0.9745 0.9749 1.0506 1.0000 1.2224 
0.0112 0.0057 0.0029 0.0014 0.0007 0.0003 
0.9745 0.9749 1.0506 1.0000 1.2224 
0.0112 0.0057 0.0029 0.0014 0.0007 0.0003 
0.9745 0.9749 1.0506 1.0000 1.2224 
0.0112 0.0057 0.0029 0.0014 0.0007 0.0003 
0.9745 0.9749 1.0506 1.0000 1.2224 
0.0112 0.0057 0.0029 0.0014 0.0007 0.0003 
0.9745 0.9749 1.0506 1.0000 1.2224 


(a0) ie) (a9) (en) 
Jl, sisi 
| 2 
iss (a 
ies) “I a oF 


T 
oO 


(en) 
lI 
SEEDED 


3 
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* Numerical Sol. 
Exact Sol. 


epsilon=1 ot 


epsilon=1 o* 


epsilon=1 of? 


Figure 1: Computational solution of the given Example 1 for the fixed value N = 100 
and various values of ¢ 


* Numerical Sol. 
0.87 Exact Sol. 


epsilon=107! 


epsilon=1 o 


epsilon=1 o3 4 


Oo 0.2 0.4 0.6 0.8 1 


Figure 2: Computational solution of the given Example 2 for the fixed value N = 100 
and various values of ¢ 
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Numerical Sol. 
0.87 Exact Sol. 


epsilon=1 ot 


epsilon=1 o? 


-0.8 wor” | epsilon=10-9 


Figure 3: Computational solution of the given Example 3 for the fixed value N = 100 
and various values of ¢ 


1.6 

- Numerical Sol. 
1.44 Exact Sol. 4d 
1.2 


epsilon=1 of 


0.8 epsilon=1 on 
0.67, 
0.47 
epsilon=10-° 
0.2 \ \ \ 1 
oO 0.2 0.4 0.6 0.8 1 


Figure 4: Computational solution of the given Example 4 for fixed value of N = 100 and 
various values of ¢ 


5 Conclusion 


We have derived an exponentially fitted finite difference tridiagonal scheme 
for solving singularly perturbed two-point BVPs at one endpoint (left or 
right). We have carried out the convergence analysis for the proposed scheme 
and performed the numerical experiments on four example problems (three 
left-end and one right-end problems) for various values of N = 1/h and 
perturbation parameter ¢, which show that the scheme is of almost first- 
order accurate. The computational results in terms of maximum absolute 
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error are presented in Tables 1—4. It is easily observed from the tables that 
the proposed method is capable of producing highly accurate results for a 
fixed value of mesh size h, when the perturbation parameter ¢€ tends to 0. 
The maximum absolute errors are becoming uniform for any fixed values of 
N when ¢ — 0. Furthermore, one can easily observe from Tables 1, 3, and 
4 that the proposed exponential fitted finite difference scheme is capable of 
producing first-order accurate uniformly convergent solution for any fixed 
value of mesh size h = 1/N when perturbation parameter € tends to 0. 
In Figures 1-4, we present our solution and the exact solution for various 
values of ¢ and a fixed value of N. Clearly, as shown in the figures, the 
numerical solution and the exact solution are very close within the boundary 
layers for smaller values of e. Notably, the novelty of our method lies in its 
independence from both deviating arguments and fitted meshes [15]. 
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